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Abstract 

In this paper the optimal control of flocking models with random in¬ 
puts is investigated from a numerical point of view. The effect of uncer¬ 
tainty in the interaction parameters is studied for a Cucker-Smale type 
model using a generalized polynomial chaos (gPC) approach. Numerical 
evidence of threshold effects in the alignment dynamic due to the ran¬ 
dom parameters is given. The use of a selective model predictive control 
permits to steer the system towards the desired state even in unstable 
regimes. 
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1 Introduction 

The aggregate motion of a multi-agent system is frequently seen in the real 
world. Common examples are represented by schools of fishes, swarms of bees 
and herds of sheep, each of them natural phenomena with important appli¬ 
cations in many fields such as biology, engineering and economy [16]. As a 
consequence, the significance of new mathematical models, for understanding 
and predicting these complex dynamics, is widely recognized. Several heuristic 
rules of flocking have been introduced as alignment, separation and cohesion 
[17, 18]. Nowadays these mathematical problems, and their constrained ver¬ 
sions, are deeply studied both from the microscopic viewpoint [5, 7, 22] as well 
as their kinetic and mean-field approximations [2, 4, 8, 9, 10]. We refer to [16] 
for a recent introduction on the subject. 

An essential step in the study of such models is represented by the intro¬ 
duction of stochastic parameters reflecting the uncertainty in the terms defining 
the interaction rules. In [1, 13, 21] the authors are concerned with the study 
of self-organized system including noise term modeling the interaction with the 
surrounding environment. In this paper we focus on the case where the un¬ 
certainty acts directly in the parameter characterizing the interaction dynamic 
between the agents. 

We present a numerical approach having roots in the numerical techniques 
for Uncertainty Quantification (UQ) and Model Predictive Control (MPC). 
Among the most popular methods for UQ, the generalized polynomial chaos 
(gPC) has recently received deepest attentions [19]. Jointly with Stochastic 
Galerkin (SG) this class of numerical methods are usually applied in physical 
and engineering problems, for which fast convergence is needed. Applications 
of gPC-Galerkin schemes to flocking dynamics, and their controlled versions, is 
almost unexplored in the actual state of art. 

We give numerical evidence of threshold effects in the alignment dynamic 
due to the random parameters. In particular the presence of a negative tail in 
the distribution of the random inputs lead to the divergence of the expected 
values for the system velocities. The use of a selective model predictive control 
permits to steer the system towards the desired state even in such unstable 
regimes. 

The rest of the article is organized as follows. In Section 2 we introduce 
briefly a Cucker-Smale dynamic with interaction function depending on stochas¬ 
tic parameters and analyze the system behavior in the case of uniform inter¬ 
actions. The gPC approach is then summarized in Section 3. Subsequently, 
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in Section 4 we consider the gPC scheme in a constrained setting and derive a 
selective model predictive approximation of the system. Next, in Section 5 we 
report several numerical experiments which illustrate the different features of 
the numerical method. Extensions of the present approach are finally discussed 
in Section 6 . 


2 Cucker-Smale dynamic with random inputs 

We introduce a Cucker-Smale type [ 8 ] differential system depending on a random 
variable 9 £ fi C K. with a given distribution f{0). Let (xi,Vi) £ R. 2d ,d > 1, 
evolving as follows 


£i(0,t) = Vi(8,t) 


N 


Vi{8,t) = A ^ ^ H ( Xi ’ X 0 ) ( v 3 (M) -Vi{0,t)) 
j =1 


( 2 . 1 ) 


where K is a time dependent random function characterizing the uncertainty in 
the interaction rates and H {-, •) is a symmetric function describing the depen¬ 
dence of the alignment dynamic from the agents positions. A classical choice 
of space dependent interaction function is related to the distance between two 
agents 


H(x,y) = 


( 2 . 2 ) 


(1 + \x - y| 2 ) 7 ’ 

where 7 > 0 is a given parameter. Mathematical results concerning the system 
behavior in the deterministic case (K = 1) can be found in [ 8 ]. In particular 
unconditional alignment emerges for 7 < 1/2. Let us observe that, even for the 
model with random inputs ( 2 . 1 ), the mean velocity of the system is conserved 
in time 


N 


N 

i—l 

since the symmetry of H implies 

N 


V(8,t) = yV(9,t) = 0, 


(2.3) 


N 


^2 H{x i ,x j )v j {e,t) = ^2 H(xi,Xj)vi{9,t). 
ij—l ij=l 


Therefore we have V(8,t) = V(0,O). 


2.1 The uniform interaction case 

To better understand the leading dynamic let us consider the simpler uniform 
interaction case when H = 1, leading to the following equation for the velocities 

Vi{0, t) = £( Wi (0, t ) - Vi {e , t)) = K(0, t)(V(8, 0) - Vi(0, t)). (2.4) 

j=i 
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The differential equation (2.4) admits an exact solution depending on the ran¬ 
dom input 9. More precisely if the initial velocities are deterministically known 
we have that 


Vi(9,t) = V+ (vi(0) - V)exp | - J K(0,s)ds}, (2.5) 

where V = V(0) is the mean velocity of the system. In what follows we analyze 
the evolution of (2.5) for different choices of K(9,t) and of the distribution of 
the random variable 9. 


Example 1 

Let us consider a random scattering rate written in terms of the following de¬ 
composition 

K(6,t) = k(9)h(t) (2.6) 

where h(t) is a nonnegative function depending on t £ R + . The expected 
velocity of the i-th agent is defined by 


Vi(t) =E g [vi(d,t)\ = / Vi(9, t)f{9)d9 
J n 

then each agent evolves its expected velocity according to 

Vi(t ) = V + (t'i(O) — V) exp | — k{9) J /i(s)ds j f(9)d9. 


(2.7) 


( 2 . 8 ) 


For example, let us chose k(9 ) = 9 , where the random variable is normally 
distributed, i.e. 9 ~ cr 2 ). Then, for each i = 1,..., N, we need to evaluate 

the following integral 


V 


'Jj (o) — v 
\f2na 2 


exp \ — 9 f h(s)ds > exp < 
i ^ Jo - 1 *• 


(g-M ) 2 

2 a 2 


\d9. 


(2.9) 


The explicit form is easily found through standard techniques and yields 

Vi{t) = V + K(0) - V) exp | - /i J h(s)ds+ a ^-^J /i(s)ds^ j. (2.10) 

From (2.10) we observe a threshold effect in the asymptotic convergence of the 
mean velocity of each agent toward V. It is immediately seen that if 

h(s)ds > ft (2-11) 

cr 2 

it follows that, for t —> +oo, the expected velocity Vi diverges. In particular, 
if h(s) e 1 we have that the solution starts to diverge as soon as t > n/cr 2 . 
Note that, this threshold effect is essentially due to the negative tail of the 
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normal distribution. In fact, if we now consider a random variable taking only 
nonnegative values, for example exponentially distributed 9 ~ Exp(A) for some 
positive parameter A > 0, from equation (2.8) we obtain 

r+oo 

€i(t) =V + ( Vi (0) - V) / e~ et X e~ x9 dd, (2.12) 

Jo 

which corresponds to 


«i(t)=V + (t>i(0)-V) 


A 

t + A’ 


(2.13) 


and therefore Vi(t) — > V as t —> oo. Then independently from the choice of 
the rate A > 0 we obtain for each agent convergences toward the average initial 
velocity of the system. Finally, in case of a uniform random variable 9 ~ U( [a, b]) 
we obtain 

«i(t) = V + K(0) - V) [ -^e~ et d9 (2.14) 

J a O CL 


that is 


Vi(t) = V + 


«i(0) - V 

b — a 




(2.15) 


which implies the divergence of the system in time as soon as a assumes negative 
values. 


Example 2 

Next we consider a random scattering rate with time-dependent distribution 
function, that is 

K(9,t) = 9(t) (2.16) 

with 9(t) ~ f(9,t). As an example we investigate the case of a normally 
distributed random parameter with given mean and time dependent variance, 
9 ~ A f(n, It is straightforward to rewrite 9 as a translation of a standard 

normal-distributed variable 9 , that is 

0 = H + cr(t)9, (2.17) 

where 9 ~ Af( 0,1). The expected velocities read 

Vi(t) = V + ^(0)_ f exp l—pt —9 ( cr(s)(is| exp < — i d9, (2.18) 

V27T ./ir *• Jo > y 2 J 

which correspond to 

Vi(t) = V + (Ui(0) - V)exp | - lit + * cr(s)dsj j. (2.19) 


5 








Similarly to the case of a time independent normal variable a threshold effect 
occurs for large times, i.e. the following condition on the variance of the distri¬ 
bution 

(/ a (s)ds^ > 2 pt (2.20) 

implies the divergence of the system (2.4). As a consequence instability can be 
avoided by assuming a variance decreasing sufficiently fast in time. The simplest 
choice is represented by a(t) = l/t a for some a £ [1/2,1). The condition (2.20) 
becomes 

/ t 1 - 01 \ 2 

( —) > 2 *. ( 2 . 21 ) 

For example, if a = 1/2 the previous condition implies that the system diverges 
for each p < 2. 

3 A gPC based numerical approach 

In this section we approximate the Cucker-Smale model with random inputs 
using a generalized polynomial chaos approach. For the sake of clarity we first 
recall some basic facts concerning gPC approximations. 

3.1 Preliminaries on gPC approximations 

Let (12, J 7 , P ) be a probability space, that is a ordered triple with any set, T 
a a —algebra and P : T —> [0,1] a probability measure on J 7 , where we define a 
random variable 

9 : (Q,T) ^ 

with £?r the Borel set of R. Moreover let us consider S C R d , d > 1 and 
[0,T] C R certain spatial and temporal subsets. For the sake of simplicity we 
focus on real valued functions depending on a single random input 

g(x,9,t) : S x n x T-* R d , g £ L 2 (fi, J 7 , P). (3.1) 

In any case it is possible to extend the set-up of the problem to a p —dimensional 
vector of random variables G = {6 1 ,..., 0 P ), see [12]. Let us consider now the 
linear space of polynomials of 6 of degree up to M, namely P m{&)- From classical 
results in approximation theory it is possible to represent the distribution of 
random functions with orthogonal polynomials {4 > fc(0)}fc/ =o , i.e. an orthogonal 
basis of L 2 (fl, J 7 , P) 


E e [$ fe (0)$ fe (0)] =E e [<f> h (6)] 2 6 hk 

with Shk the Kronecker delta function. Assuming that the probability law 
P(9~ 1 (B)),\/B £ Br, involved in the definition of the introduced function 
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Table 1: The different gPC choices for the polynomial expansions 


Probability law of 9 

Expansion polynomials 

Support 

Gaussian 

Hermit e 

(- 00 , +oo) 

Uniform 

Legendre 

[a, b] 

Beta 

Jacobi 

[a, b] 

Gamma 

Laguerre 

[0, Too) 

Poisson 

Charlier 

N 


g{x 1 6 , t) has finite second order moment, then the complete polynomial chaos 
expansion of g is given by 

9(x,0,t) = ^ g m (x,t)$ m (9). (3.2) 

m£ N 

Accordingly to the Askey-scheme, result that pave a connection between ran¬ 
dom variables and orthogonal polynomials [19, 20], we chose a set of polynomials 
which constitutes the optimal basis with respect to the distribution of the in¬ 
troduced random variable in agreement with Table 1. 

Let us consider now a general formulation for a randomly perturbed problem 

V{x,t,6;g) = f(x,t,9) (3.3) 

where we indicated with T> a differential operator. In general the randomness 
introduced in the problem by 8 acts as a perturbation of T >, of the function g or 
occurs as uncertainty of the initial conditions. In this work we focus on the first 
two aspects assuming that initial positions and velocities are deterministically 
known. 

The generalized polynomial chaos method approximate the solution g(x, 9 , t ) 
of (3.3) with its Mth order polynomial chaos expansion and considers the 
Galerkin projections of the introduced differential problem, that is 

E e [V(x,t,9-,g)-<f> h (9)\=E e \f(x,9,t)-$ h (9)\, h = 0,1,..., M. (3.4) 

Due to the Galerkin orthogonality between the linear space Pm and the error 
produced in the representation of g(x,9,t) with a truncated series, it follows 
that from (3.4) we obtain a set of M + 1 purely deterministic equations for 
the expansion coefficients g m (x,t). These subproblems can be solved through 
classical discretization techniques. From the numerical point of view through a 
gPC-type method it is possible to achieve an exponential order of convergence 
to the exact solution of the problem, unlike Monte Carlo techniques for which 
the order is 0(l/yM) where M is the number of samples. 

3.2 Approximation gPC of the alignment model 

We apply the described gPC decomposition to the solution of the non-homogeneous 
differential equation Vi(9, t) in (2.5) and to the stochastic scattering rate K[9 , t), 
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(3.5) 


i.e. 


M 


M 




(M) = Vi, m (t)$m(0), K M (0,t) = ££,(t)$ z ( 0 ) 


m =0 


1=0 


where 

Vi,m{t) = E e [vi{6,t)$ m (6)} ki(t) = E e [K(0,t)®i{9 )], (3.6) 

we obtain the following polynomial chaos expansion 
d M i N m 

22 Vi,m®m(0) =^^2 H {Xi,Xj) 22 K l{t){Vj, m - Uj, m )$Z (0)$ m (0)- 


dt 


m =0 


j=i 


l,m=0 


Multiplying the above expression by an orthogonal element of the basis <&h(0) 
and integrating with respect to the distribution of 9 


Eg 


M 


22 


m —0 


= Efl 


AT 


M 


y: Ki(t){vj,rn-v i , rn )$i(6)$ m (o)$ h (6) 

j =1 l,m =0 

we find an explicit system of ODEs 

JV M M 

- , \ 1 
7).- u(t) = - 

dt 


^V iih {t) = N ^ || 2Ar 22 H ( x i’ x i) 5Z ~ W™) 51 Ki(t)eimh 


\<S>h\\ 2 N 


3=1 


m —0 


Z=0 


1 


AT 


M 


(3.7) 


— ^y- 5 ) Xj') 5 . (Jtj.m 

j—1 m =0 


where e Jm?l = E e [$j(0)$ m (0)$ h (0)] and 


M 


k mh (t) ~ 19 5 1 ki(t)eimh- 


;=o 


We recall that the gPC numerical approach preserves the mean velocity of the 
alignment model (2.4). In fact, from (3.7) follows 


M 


N d i N 

22 5Z H ( X i’ X i) 51 Krnh{t)Vj,r, 


dt ' K ' N 

2—1 j,2—1 m—0 


AT 


M 


(3.8) 


^ ^ H(Xi,Xj) ^ ^ ^rnh^)^i,rn — b? 
j,2=l ra=0 

thanks to the symmetry of i/. More generally it can be shown that if 

AT 


N 


22 v i(8,t) = v> 


2=1 








where V is time-independent, then its gPC decomposition is also mean-preserving 
since 

N M M 

jv £ £ E « k(M)<M0)] <M0) = £ 

2=1 m=0 m —0 

M 

= V^E e [l-$ m (0)] $ m (0)=V. 

m—0 

Remark 1. TTie gPC approximation (3.7) can be derived equivalently without 
expanding the kernel function K(9,t). In this way one obtains 

d 1 N M 

^ ' (Vj,m ^i,rn)I^Tnh (3.9) 

j—1 m—0 

where now 

K mh (t) = — l -^E 0 [K(O,t)$ m <f> h }. 

Note that, since in general N M, the overall computational cost is 0(MN 2 ). 


N 


— ^2vi{6,t)$m{0) 


®m{0) 


4 Selective control of the gPC approximation 


In order stabilize the gPC approximation of the Cucker-Smale type model (2.1) 
with random inputs, we introduce an additional term which acts as control of the 
approximated dynamic. More specifically we modify the approximation of the 
alignment model ( 2 . 1 ) by introducing a control term Uh to the ht h component 
of its gPC approximation 

d 1 N M 

(0 =-^y ^H(Xj,Xj) jT Kmh(t)(Vj,m(t ) - Vi,m(t)) + U h Q{Vi,h) (4.1) 

j—1 m—0 


where Uh is a solution of 


Uh = arg min 



Vd,h) 2 dt 


V 

2 



u h {t) 2 dt 


(4.2) 


where v > 0 is a regularization parameter and {vd,o, Vd,i, ■ ■ ■ j v<i,m) are the 
desired values for the gPC coefficients. For example 


Vd,h = ^e[vd®h(0)] = ^[ 4 ^( 60 ] 


v d h = 0 
0 h = 1,..., M, 


(4.3) 


where Vd is a desired velocity. 

Moreover the controller action is weighted by a function, Q(-), such that 
Q{vi,h) £ [— L,L],L > 0. Due to the dependence of the controller effect from 
the single agent velocity, we refer to this approach as selective control, see [3]. 


9 







In order to tackle numerically the above problem, whose direct solution is 
prohibitively expansive for large numbers of individuals, we make use of model 
predictive control (MPC) techniques, also referred to as receding horizon strat¬ 
egy or instantaneous control [15]. These techniques has been used in [2, 3, 4] in 
the case of deterministic alignment systems. 


4.1 Selective model predictive control 


The basic idea is to consider a piecewise constant control, 

m— 1 

Uh(t) = 

n —0 

on a suitable time discretization. In this way is possible to determine the value 
of the control ttJJ G M, solving for a state Vi t h the (reduced) optimization problem 


d 1 1V1 

= ^ ^ ^ H(Xi, Xj ) ^ ^ “1“ 

j—1 m —0 

Vi,h(t n ) = Vi,h , 

l>t n+1 l N v \ 

K = arg min I — ^ - v d ,h) 2 + ^ u h (t) 2 \ ds,u% G [uh,L,K,Rl 

(4.4) 

Given the control on the time interval [t n ,t n+1 ], we let evolve Djj, according 
to the dynamics 


d i JV M 

= ~jy ^ ^ ( -Kmh (t) (hj',i 


i W - Vi,m(t)) + KQ(v* h (t)) (4.5) 


i=i 


m=0 


in order to obtain the new state ^ = fi ? ;y i (t" +1 )■ We again solve (4.4) to obtain 
u£ +1 with the modified initial data and we repeat this procedure until we reach 
nAt = T. 

The reduced optimization problem implies a reduction of the complexity 
of the initial problem since to an optimization problem in a single real-valued 
variable u%. On the other hand the price to pay is that in general the solution 
of the problem is suboptimal respect to the full one (4.1)-(4.2). 

The quadratic cost and a suitable discretization of (4.5) allows an explicit 
representation of v% in terms of Vi t h and , as a feedback controlled system 
as follows 


At 


N 


M 


Kt 1 = + jy E H l E %kh(v?,m ~ Km) + A tulQl 


hi 


j—1 m =0 


Vi,h = V i,h, 

At 


N 




i —1 


(4.6) 
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where //”■ = H(x™, x ") and Q™ h = Q(v™ h ). Note that since the feedback control 
uj) in (4.6) depends on the velocities at time n+ 1, the constrained interaction at 
time n is implicitly defined. The feedback controlled system in the discretized 
form results 


u n + 1 

i,h 


3 =1 


N 


M 

m =0 


At 2 


N 

I'- vim) - ^ E 1 

i=i 


vU 1 V d ,h)Qj,hQi,hi 


v i,h = v i 


h • 


Again the action of the control is substituted by an implicit term representing 
the relaxation toward the desired component of the velocity Vd,h, and it can be 
inverted in a fully explicit system. 

Considering the scaling for the regularization parameter v = fcAt, the pre¬ 
vious scheme is a consistent discretization of the following continuos system 


d 

dt 




1 N M 

j =1 m—0 


1 x A 

+ -r; “ v j!h (t))Q(v jjh (t))Q(vi, h (t)). 


(4.7) 


Now the control is explicitly embedded in the dynamic of the hth component of 
the gPC approximation as a feedback term, and the parameter k > 0 determines 
its strength. 


4.2 Choice of the selective control 

For the specific choice of weight function Q{-) = 1 we refer in general to non 
selective control. Note that in this case the action of the control is not strong 
enough in order to control the velocity of each agent, indeed in this case we are 
able only to control the mean velocity of the system. In fact the control term 
is reduced to 

-(t>d,h-V/.), (4.8) 

K 

where Vh is the h —th coefficient of the expansion of V, that is 

1 N 

Vh = 

3 =1 

Then, only the projections of the mean velocity are steered toward the respective 
components of the target velocity, i.e. as soon as k —> 0 it follows that Vh = Vd,h- 
Therefore, the choice of the selective function Q(-) is of paramount importance 
to ensure the action of the control on the single agent. 
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In principle one can address directly the control problem on the original 
system (2.1) as 


Xi{0,t) = Vi(9,t ) 

I Vi(9,t) = K ^ ^ Y^H(x h xj) (vj (9, t) - Vi(9,t)) + uQ(vi(9,t)), 
where the control term u is solution of 


(4.9) 


u = arg min 

u 


I Jj'52( v i( 9 ’ t )- v d) 2dt + 7;[ u(t) 2 dt 

Jo v i=1 Jo 


(4.10) 


Here Vd G M. d is a target velocity and v > 0 a regularization parameter. Similarly 
to previous subsection, through the approach presented in [2, 3, 4], we can derive 
the time-continuos MPC formulation which explicitly embed the control term 
in the dynamic as follows 


±i(9,t)= Vi(9,t) 

Vi{9,t) = ^^H( x i, x j)( v j(9,t) — Vi(9,t)) 

j =i 

1 N 

_ v j( 6 ' i )) ( 5(^( i > e ))Qi v i{^ 6 ))- 

3 = 1 


(4.11) 


Now the gPC approximation of (4.11) can be obtained as in Section 3 and leads 
to the set of ODEs 


d 

dt 




1 N M 

= -jy ^ > H{Xi, Xj ) ^ ^ 
j—1 m =0 



^R h {vf ,vf), 


3= 1 


(4.12) 


where 

Rh(v™,vY) = [(vd - vf)Q (v™(0,t)) Q (uf (M)) <M0)] . (4.13) 


In general systems (4.12) and (4.7), without further assumptions on the selective 
function Q(-), are not equivalent. In addition to the non selective case, there 
exist at least one choice of selective control that makes the two approaches 
totally interchangeable. In fact, taking 


QM = 7 * 1 =, (4.14) 

v££f= 


12 












Figure 1: Sketch of the two numerical approaches to solve the control problem 
with uncertainty, combining MPC and gPC. In both cases, of non selective 
control, i.e. Q(-) = 1, and of selective control with Q(-) defined in (4.14) the 
two approaches are equivalent. 


we have Q(-) £ [—1,1] and the control term in (4.12) takes the following form 

1^1 1 

= K ||$ ||2 Ee = - (v d , h - v ilh ). (4.15) 

j =i 

Similarly the control term in (4.7) reduces to 

1 N 1 
^ - Vj, h (t))Q(vj,h{t))Q(vi,h(t)) = - (y d ,h ~ Vi,h) , (4.16) 

j =i 

and therefore system (4.12) coincides with (4.7). Note that as k —> 0 both 
systems are driven towards the controlled state Vi,h = Vd,h which implies a 
strong control over each single agent. 

In Figure 1 we summarize the two equivalent approaches. In the case of non 
selective control and of selective function given by (4.14) the constrained gPC 
system can be obtained from our initial unconstrained model (2.1) through two 
different but equivalent methods. The first approximates the solution of the 
Cucker-Smale type model via the gPC projection and then introduces a con¬ 
trol on the coefficients of the decomposition through a MPC approach in order 
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to steer each component to (vd,o, «d,i, • • •, VcI,m)- Whereas the second method 
considers a constrained Cucker-Smale problem (4.9), introduces its continuous 
MPC approximation and then computes the gPC expansion of the resulting 
system of constrained differential equations. 

Remark 2. We remark that the choice of Q(-) stated in (4.14), for which the 
two approaches sketched in Figure 1 are identical, is equivalent to consider the 
constrained dynamic (4.9), modified as follows 


±i(0,t ) = Vi(9,t) 

Vi(0,t) = YHjxi'XMvMt) - Vi(0,t )) + Ui, 

j=i 


(4.17) 


where the control term, Ui for each agent i = 1 is given by the mini¬ 

mization of the following functional 


\ (' I \ ^ ' - 

J(v 1 ,...,v N -,ui,...,u N ) = - / — (vi(9,t) 

J 0 ; _ 1 


~ V d ) 2 + 


dt. (4.18) 


Since the functional is strictly convex, applying the (MPC) procedure to (4.17)- 
(4.18), we obtain as first order approximation for the solution of the optimal 
control problem the feedback control term 


Ui = ~(vd - Vi), i = 1,... ,1V. 

K 


(4.19) 


Thus the same considerations on the equivalence of the approaches hold. 


5 Numerical tests 

We present some numerical experiments of the behavior of the flocking model 
in the case of a Hermite polynomial chaos expansion. This choice corresponds 
to the assumption of a normal distribution for the stochastic parameter in the 
Cucker-Smale type equation (2.1) and in its constrained behavior (4.7). Nu¬ 
merical results show that the introduced selective control with the weight func¬ 
tion (4.14) is capable to drive the velocity to a desired state even in case of 
a dynamic dependent by a normally distributed random input, with fixed or 
time-dependent variance. In the uniform interaction case, since the effect of 
agents’ positions do not influence the alignment we report only the results of 
the agents’ velocities. 

5.1 Unconstrained case 

In Figures 2 and 3 we present numerical results for the convergence of the error 
using the gPC scheme described in equation (3.7) for H = 1 and solved through 
a 4th order Runge-Kutta method. In particular Figure 2 shows the behavior 
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Figure 2: Error convergence for increasing number of polynomials in the gPC 
decomposition approximation. Left: convergence of the mean error at two fixed 
times T = 1 and T = 5. Right: convergence of the variance error. In both 
cases we considered a random time-independent scattering K(9,t) = 9, where 
the random variable 9 is normally distributed 7V(2,1/2). The system of ODEs 
is solved through a 4th order Runge-Kutta with At = 10~ 5 . 


of the error with respect to increasing terms of the gPC decomposition. Here 
we considered the average in time of the error for the mean and the variance at 
time t > 0 in the L 1 norm 


Ey{t) 


if 

N 


N 


Es 2 W ~ tv 


i—1 






(5.1) 


where 


af(t) = Eg - Vi(t)Y 


(5.2) 


with i>i(9,t) and Vi(t) defined in (2.5) and (2.7). Observe that if the scattering 
rate K{9 1 t) is of the from described in (2.6) with h(-) = 1 and k{9) ~ Af(fi,a 2 ) 
than, in addition to the explicit evolution for the expected velocity as in (2.8), 
we can obtain the exact version for the evolution of the variance of the ith agent 


of (f) = (uj(0) — V) 2 (exp{— 2fit + 2 a 2 t 2 } — exp{— 2 fit + cr 2 t 2 }) . (5.3) 


In (5.1) we indicated with 


a i’ M (t) the approximated variance 


M 

= Y, vlh(t)^e[$h(0) 2 ] - vf, 0 (t). (5.4) 

h =0 

It is easily seen how the error decays spectrally for increasing value of M, how¬ 
ever the method is not capable to go above a certain accuracy and therefore 
for large M a threshold effect is observed. This can be explained by the large 
integration interval we have considered in the numerical computation, and by 
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Figure 3: Evolution of the variance-error E v {t) defined in equation (5.1) for 
the gPC decomposition for the unconstrained model (2.4) with K(9,t) = 9 ~ 
A/"(2,1/2) over the time interval [0,T] with T = 5 and time step At = 10 -5 . 


the well-known loss of accuracy of gPC for large times [12]. In the case of the 
error of the variance, Figure 3, the gPC approximation exhibits a slower con¬ 
vergence with respect to the convergence of the mean. Next in Figure 4 we see 
how for large times the solution of the differential equation (2.4) diverges and 
the numerical approximation is capable to describe accurately its behavior only 
through an increasing number of Hermite polynomials. 


5.2 Constrained uniform interaction case 


In Figure 5 we show different scenarios for the uniform interaction dynamic 
with constraints. In the first row we represents the solution for N = 10 agents, 
whose dynamic is described by equation (4.7) with Vd = 1, different values of k 
originate different controls on the average of the system, which however do not 
prevent the system to diverge. In the second row we show the action of selective 
control (4.14). It is evident that, with this choice, we are able to control the 
system also in the case with higher variance. 

Observe that the numerical results are coherent with the explicit solution 
of the controlled equation. Let us consider the time-independent scattering 
rateA'(0, t) = 6 ~ A f(fi, a 2 ), then from the equation 

4; Vi(0,t ) = 0(V - Vi(8,t)) + ~(v d - Vi(0,t)) (5.5) 

at k 

we can compute the exact solution given Vj(9. 0) = Uj(0) 


Vi(9,t) 


kV 9 + Vd 
k9 + 1 



kV 9 + Vd 
k9 + 1 


ex P i - 


{~( 0+ £W 


(5.6) 


The asymptotic behavior of the expected value of (5.6) can be studied similarly 
to what we did in Section 2.1. In other words in order to prevent the divergence 
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Figure 4: Left: 6th order Hermite gPC decomposition solved through a 4th order 
Runge-Kutta. Right: 10th order Hermite gPC decomposition solved through a 
4th order Runge-Kutta. In both cases the final time considered is T = 6, with 
time step At = 10 -2 . 


of the leading term of the controlled expected exact solution we might study 

exp {_( /1+ I) t+ ^} j (5.7) 

which diverge if 

+ (5 - 8) 

Then for each fixed time we could select a regularization parameter n > 0 so 
as to avoid the divergence of (5.6). Moreover we can observe that in the limit 
tv —^ 0 the introduced selective control is capable to correctly drive the system 
(5.5) for each t > 0. 

In Figure 6 we consider the system with random time-dependent scattering 
rate 9 ~ J\f(n,cr 2 (t)). The dynamic shows how, for the choice of time depen¬ 
dent variance described in Remark 2.1, that is cr(t) = 1 /s a with a = 1/2, the 
convergence depends from the mean value of the random input. In particular 
numerical experiments highlight the threshold effect for /i = 2 which we derived 
in Section 2. In the second figure we show that the action of the selective control 
(4.14), with desired velocity vj = V, is capable to stabilize the system and drive 
the velocities towards the desired state. 

5.3 Constrained space dependent case 

Next let us consider the full space non homogeneous constrained problem (2.1) 
with the interaction function defined in (2.2). In this case we assume that 
K{9) = 6 with 9 ~ 7V"(//, cr 2 ). In Figure 7 and 8 we consider a system of N = 100 
agents with Gaussian initial position with zero mean and with variance 2 and 
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Gaussian initial velocities clustered around ±5 with mean 1/10 . The numerical 
results for (3.7) have been performed through a 10th order gPC expansion. The 
dynamic has been observed in the time interval [0, 5] with At = 10 -2 . In Figure 
8 we see how the selective control is capable to drive the velocity of each agent 
to the desired state V&. In fact in case of no control, see Figure 7, we have that 
the velocities of the system naturally diverges. 

6 Conclusions 

We proposed a general approach for the numerical approximation of flocking 
models with random inputs through gPC. The method is constructed in two 
steps. First the random Cucker-Smale system is solved by gPC. The presence of 
uncertainty in the interaction terms, which is a natural assumption in this kind 
of problems, leads to threshold effects in the asymptotic behavior of the system. 
Next a constrained gPC approximation is introduced and approximated though 
a selective model predictive control strategy. Relations under which the intro¬ 
duction of the gPC approximation and the model predictive control commute 
are also derived. The numerical examples illustrates that the assumption of pos¬ 
itivity of the mean value of the random input is not sufficient for the alignment 
of the system but that a suitable choice of the selective control is capable to sta¬ 
bilize the system towards the desired state. Extension of this technique to the 
case of a large number of interacting agents through mean-field and Boltzmann 
approximations are actually under study. 
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Evolution with 9 ~ Af(2,1) 



Evolution with 9 ~ M(2, 0.5) 



Evolution with 9 ~ A/”(2,1) Evolution with 9 ~ N(2, 0.5) 



Figure 5: Evolution of the uniform interaction alignment model (4.7) with N = 
10 agents, at t = 0 distributed around V = 2 with unitary variance, depending 
on a normal random parameter. Left column: 8 ~ 7V"(2,1). Right column 
9 ~ 7^(2, 0.5). The control term shows its ability to steer the system towards 
desired velocity Vd = 1, with different intensities k = 1 and k = 0.1, when 
k = oo the control has no influence. First row shows the action of the control 
acting just on the average velocity, Q = 1. Second row shows the action of 
selective control with Q(-) as in (4.14). 
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Figure 6: Solution of the uniform interaction case with time dependent random 
parameter 0 distributed accordingly to a normal distribution A/"(/x, cr 2 (t)), with 
a time-dependent standard deviation a(t) = l/f Q , and a = 1/2. Left: we see 
the threshold for different values of /i, i.e. for ^ < 2 the system diverges. Right: 
solution of the constrained model with k = 0.1, observe that we are able to steer 
the system to the desired velocity Vd = V, i.e. the initial mean velocity of the 
system, using the selective control described in (4.14). 
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Position 
(a) t=0 


(c) t=2 


Position 

(b) t=l 


(d) t=3 



Position 


(e) t=4 



Figure 7: Numerical solution of (4.7), with 7 = 0.05 < 1/2, ( = 0.01, through 
a 10th order gPC Hermite decomposition (3.7) with k = 00 with time step 
At = 10~ 2 . The random input is normally distributed 6 ~ A/”(2,1). 
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Position 

(a) t=0 


(b) t=l 



Position Position 


(c) t=2 


(d) t=3 


(e) t=4 


(f) t=5 


Figure 8: Numerical solution of (4.7), with 7 = 0.05 < 1/2, £ = 0.01, through a 
6 t/i-order gPC Hermite decomposition for the selective control (3.7) with time 
step At = 10 -2 . Here we considered a normally distributed random input 
0 ~ Af{2, 1), the desired velocity is = 0 and the control parameter is n = 1. 
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